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Networks of biofilaments are essential for the formation of cellular structures that support various 
biological functions. For the most part, previous studies have investigated the collective dynamics 
of rod-like biofilaments; however, the shapes of the actual subcellular components are often more 
elaborate. In this study, we considered an active object composed of two active filaments, which 
represents the progression from rod-like biofilaments to complex-shaped biofilaments. Specifically, 
we numerically assessed the collective behaviors of these active objects in two dimensions and ob¬ 
served several types of dynamics depending on the density and the angle of the two filaments as 
shape parameters of the object. Among the observed collective dynamics, a moving density band 
that we named a ‘moving smectic’ is introduced here for the first time. By analyzing the trajectories 
of individual objects and the interactions among them, this study demonstrated how interactions 
among active biofilaments with complex shapes could produce collective dynamics in a non-trivial 
manner. 

PACS numbers: 87.16.Ka 05.65.+b 


I. INTRODUCTION 


Organized networks of biofilaments such as micro¬ 
tubules and actin fibers play essential roles in the for¬ 
mation, maintenance, and alteration of cell structures 
[ l|. For example, cortical microtubules in plant cells ex¬ 
hibit aligned structures called ‘bundles’, which mechani¬ 
cally support the cells These subcellular structures 

composed of biofilaments are usually formed through ac¬ 
tive processes driven by energy use in hydrolysis. Poly¬ 
merization and depolymerization cause length changes 
to the microtubules and cause them to collide with each 
other, which ultimately leads to their global alignment. 
Another active mechanism involves interactions between 
biofilaments and molecular motors. When the molecu¬ 
lar motors bind to the filaments and march along them, 
active forces act on the filaments and drive the forma¬ 
tion of spatio-temporal filament patterns IM3- To date, 
a number of in vitro experiments have revealed various 
types of structures resulting from this motor activity P- 
Il3l |. For example, filaments form locally ordered patterns 
such as ray-like asters and vortices, depending on the 
type and concentration of molecular motors present [Ti¬ 
ll 7]. In addition, the collective motion of biofilaments 
can emerge via their active interaction with each other, 
which has led to the interesting paradigm of ‘active mat¬ 
ter’ dm. 

Many theoretical studies have attempted to address 
the self-organization of active filaments. Based on the 
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contracting force of active filaments driven by molecular 
motors, coarse-grained and continuum kinetic equations 
were proposed to explain the inhomogeneous accumula¬ 
tion of filaments (i.e., bundles) 24, 25[ and local struc¬ 
tures such as vortices and asters [26[ ■ Similar global pat¬ 
terns were observed in another continuum model in which 
nematic collisions were taken into account [27l - l3ll |. To 
explain the emergence of vortices and asters from micro¬ 
scopic processes, Aranson and Tsimring introduced sim¬ 
ple stochastic rules for inelastic collisions of biofilaments, 
and thus derived continuum equations by which trans¬ 
port coefficients could be related to microscopic physical 
quantities [32, 33] . Aranson’s group performed Monte 
Carlo simulations of the elementary alignment processes 
to verify their theory (34| ; however, a correspondence 
with continuum theory remains elusive. In other stud¬ 
ies simulating the dynamics of filaments, new patterns, 
including stripes, were found (35|, and the viscoelasticity 
of the networks of biofilaments was discussed [3(|. 

The majority of previous studies mainly considered ob¬ 
jects with rod-like shapes as the simplest biofilaments. 
However,in vivo, filaments often bind to each other with 
additional related molecules and form molecular com¬ 
plexes, which constitute a functional unit in various bio¬ 
logical processes. For instance, microtubules on mitotic 
spindles elongate radially from a pair of centrosomes, 
which separate a mother cell into two daughter cells ]37j . 
Beneath the apical membrane of a multi-ciliated cell, the 
root of the cilium, i.e., the basal body, has an appendage 
known as the basal foot that functions as a microtubule 
organizing center. Microtubules are generated from these 
basal feet and subsequently connect with each other to 
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become organized as a cell-sized network, which is sug¬ 
gested to be involved in the ordered alignment and di¬ 
rection of the beating cilia dSj!, [39j. In this example, 
microtubules that are pivoted by the basal body can act 
as a functional unit, and their interaction may lead to 
alignment of the cilia. However, the mechanisms under¬ 
lying the function of such biofilament complexes, espe¬ 
cially the relationship between the shape of the molecular 
complexes and their emerging dynamics, remain largely 
unknown. In addition, current developments in nanobio¬ 
engineering, such as optical manipulation [40], have en¬ 
abled the design of biofilament complexes that can give 
rise to new types of the self-organization of biofilaments. 
Such technological progress has also motivated us to ex¬ 
plore the possible dynamics of biofilaments when they do 
not have a simple rod-likc shape. 

To gain fundamental insights into these mechanisms, in 
the present study, we performed numerical simulations of 
molecular complexes composed of two filaments, named 
‘active filament complexes (AFCs)’. This model rep¬ 
resents a simple investigative progression from rod-like 
filaments to complex-shaped filaments. In our model, 
the active interactions of filaments via molecular mo¬ 
tors are considered in two dimensions, whereas collisions 
and excluded volume interactions are ignored. We dis¬ 
covered interesting dynamics depending on the density 
and shape of AFCs, including a newly identified dynamic 
form termed a ‘moving smectic’ (described in detail later) 
that has not been previously reported. By tracking tra¬ 
jectories and investigating AFC interactions, we also ex¬ 
plored the mechanisms by which the various dynamic 
patterns observed might arise in the system. 

The remaining sections of this paper are organized as 
follows. The detailed model is described in Sec. HU In 
Sec. EH the observed dynamics in numerical simulations 
are classified according to quantities such as ferromag¬ 
netic (polar) and nematic (apolar) order parameters. By 
tracking individual AFCs and characterizing the cross¬ 
ing of filaments, analyses are described from Sec. IIVI to 
Sec. lVIl Finally, a comparison with earlier studies and the 
potential for future experiments are discussed in Sec. lVIll 


II. MODEL 

A. Dynamics of single active filaments 

To model the dynamics of sing le active filaments, the 
formulation given by Tanase [4J| was modified. In this 
formulation, a biofilament is described as a single rod 
diffusing in a viscous liquid. Each filament possesses po¬ 
larity associated with its elongating direction from the 
minus to plus end. Therefore, 3 N variables are sufficient 
to represent the configuration of the microtubules, i.e., 
their center position, = ( Xi,yi ), and the direction of 
their plus end, V’i (G [0, 27 t]) [Fig. Ufa)]. 

In a viscous environment where hydrodynamic interac¬ 
tions are negligible, the dynamics of a rod-like filament 


are governed by R, = C 1 F i and = ( r „ t [xF iy - yF ix ). 
Here, = (FA,, Fi y ) is the external point force acting on 
the filament, and x = ( x , y ) is the relative vector point¬ 
ing from the center of mass of the filament to the point 
of action. £ _1 = Cj^UjUi + CJ 1 (I — UjUj) represents 
the reciprocal matrix of orientation-dependent drag coef¬ 
ficients, where Uj = (cos0j, sin(9i) is the unit orientation 
vector parallel to the filament. £||, £j_, and C, ro t are drag 
coefficients associated with the longitudinal, transversal, 
and rotational motions of filaments, respectively. These 
coefficients in dilute systems are related as £|| = £j_/2, 

£rot = yrC||; where £ is the filament length and constant 
c = 6 is independent of £ at the limit of the high aspect 
ratio of the rods [4^, S^[. We fixed t by ignoring polymer¬ 
ization/depolymerization and approximate F = (A ~ (x, 
since their difference is only a factor of 2 [33 • Then, the 
dynamics of R^ has the simple form R, = £ -1 Fi. 

The point force to the i-th filament, F,;, derives from 
two independent parts. One is the thermal noise trig¬ 
gered by the random motion of the solvent molecules. 
The other is the motor-mediated interaction with other 
filaments, which is our main point of interest, formulated 
as follows. In the cell cytoplasm or solvents, many molec¬ 
ular motors such as kinesin and dynein bind to the fila¬ 
ments. While kinesin binds, a sliding motion occurs along 
the filament from the minus to the plus end. Because 
most molecular motors bind to the filaments by multiple 
binding heads, when two filaments are in close proximity 



FIG. 1. Schematic representations of the interaction between 
filaments and active filament complexes (AFCs). (a) Con¬ 
figuration of crossing single filaments. Red points represent 
the center of filaments, (b) Possible alignment of crossed fila¬ 
ments (‘zippering’). (c) AFC with shape angle 9 a and length 
£ nm. Each filament is named as an L- or R- filament. The 
center of mass in the AFC is located on the bisecting arrow 
(red). The cyan point denotes the center of mass of an AFC 
in the light-anchor model, (d) Interaction and possible align¬ 
ment of crossed AFCs. 
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and are simultaneously bound to the same motors, the 
siblings of the motors result in mechanical torque that 
orientates the two filaments by generating tension T m 
at the bindin g po int R c , a process known as ‘zippering’ 
[Fig- mb)] l44l ]. In this process, molecular motors dif¬ 
fuse ^100 times faster than biofilaments, and can hence 
be regarded as uniformly distributed [33[. In addition, 
under the high concentration of molecular motors, the 
binding to filaments and sliding along them occur at a 
constant rate. These conditions enable us to treat T m as 
a fixed parameter without explicitly considering motor 
dynamics. We assume that the direction of the zippering 
force generated by the molecular motors only depends 
on the relative angles between a contacting pair of two 
filaments [Fig. [TJa) and (b)] [4l|. Taken together, for a 
two-dimensional system in which crossings of filaments 
are not prevented (see below), the equations of motion 
for the *-th filament are expressed as 
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The first terms of the right-hand side of the two equa¬ 
tions represent ‘zippering’, and the second terms rep¬ 
resent thermal noise. Cj represents a set of filaments 
that interact with the i-th filament via motors. If 
none of the filaments interacts, the first terms van¬ 
ish. r/i(t) is a state-dependent random force obeying 
r?i(t) = 77] 1 (t)Ui + rj±(t)Vi, where V, = (-sin0j,cos0j) 
is the unit vector perpendicular to the *-th filament. 
7/|, Vi~-> and Vl ot represent Gaussian white noise follow¬ 
ing the statistics: (rjf (t)) = (■ rj±(t )) = (?y[ ot (t)) = 0, 

(v! {t)v! 00) = 2D^S{t-t'), (vH^vH*')) = 2D ± S(t-t’), 

( 1 li(t)Vi 0t {t')) = 2D rot 6(t - t'). D ||, D±, and D rot are 
longitudinal, transversal, and rotational diffusion con¬ 
stants, and are related to drag coefficients by the Ein¬ 
stein relations as D = D || = D± = ksTC, -1 and 
Drot = ksTCroti respectively, where ksT should be in¬ 
terpreted as the effective temperature. In the equations, 
total momentum is conserved if thermal noise is elimi¬ 
nated. 

In the derivation of the equations above, we ignored 
the excluded volume effect for simplicity, although it may 
contribute to apolar alignment [291, [3l[. This can be jus¬ 
tified with two assumptions [Fig. [TJb)]. One is that we 


set a pseudo two-dimensional system whose length along 
the z-axis is around 0.1 —1.0 /im. Thus, filaments cannot 
lean toward the z-direction due to their length (> 1 fj ,m 
for a microtubule), and their crossings are not prevented 
by physical contacts along the x — y plane because both 
the diameter of the filaments (~ 26 nm for a microtubule) 
and the size of the molecular motors(~80 nm for kinesin) 
are sufficiently small. The other assumption concerns the 
motor property. If the density of the motors and their 
affinity to the filaments are both high, the filaments are 
always associated with multiple motors that cross-link 
and drive the zippering interaction. Thus, in our setup 
of two-dimensional simulations, motor-mediated tension 
always acts at the crossing point and zippering will have 
a greater effect than collision. This situation holds for a 
microtubule network beneath the apical membrane and 
is consistent with many experimental setups. 


B. Model for AFCs 


To study the dynamics of active filaments with shape, 
the model equations for active single filaments were ex¬ 
tended as follows. As a simple AFC, we consider an ob¬ 
ject to which two filaments are anchored on their minus 
ends at an angle of 6 a [Fig. [He) blue point]. The lengths 
of the two filaments are set as t. Hereafter, 9 a is referred 
to as the shape angle, and the left and right filaments 
are called L- and R- filaments, respectively. We assume 
that the mass of the anchoring object is negligibly small 
compared to that of the filaments, and the motion of the 
two filaments dominates the dynamics of AFCs. For an 
AFC composed of the *-th and fc-th filaments, the cen¬ 
ter of mass of the AFC becomes R m = (Rj + Rfc)/2. 
When an AFC interacts with other AFCs, the equations 
of motion for the center of mass and the angle, R m and 
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ip m =ipi = ipk, become [Fig. QJd)], 
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comparable to those used in previous studies [U, . At 

these values, the noise is quite weak, and thus the dy¬ 
namics of the systems are mostly dominated by active 
forces. The simulation time is taken as t = 1000 s, which 
is sufficiently long compared to the relaxation times of 
all order parameters (< 100 s) introduced in Sec. IIIII 
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where ( k i) represents the counterpart term to the 
first half term of each equation, permuating the indices 
i, k, and Ci into k, i, and Ck , respectively. Here, C,; and 
Ck are sets of filaments that interact with the i-th and 
fc-th filaments, respectively. 

This model is referred to as the ‘light-anchor model’. 
In general, when the mass of filaments is mj and that 
of the anchor is m a , the center of mass of the AFC is 
positioned at m E+ 2 mE distance from the anchor along 
the bisection of the two filaments [Fig. [He) red line]. 
The light-anchor model is justified in the case of m a 
2to/. If noise is absent, the model also conserves total 
momentum and is classified into an ‘active wet system’ 
according to Ref. Q. 

Numerical simulations of AFCs were performed in two- 
dimensional systems with periodic boundary conditions. 
The initial states of AFCs are set to obey spatially and 
orientationally uniform distributions. The system size is 
fixed at L x L = 4.0 x 10 2 pm 2 unless mentioned oth¬ 
erwise. To evaluate the effect of the AFCs’ shape on 
the collective dynamics, the shape angle 6 a is controlled. 
At 9 a = 0°, the AFC coincides with a polar filament, 
whereas at 9 a = 180°, the AFC is a bipolar filament with 
length 2 i ^m. The density of the AFCs is also controlled 
by changing the number of AFCs, N. The other param¬ 
eters are fixed as follows: filament length, t = 2.0 /rm; 
translational drag coefficient, ( = 0.4x 10 -6 kg • s -1 ; ef¬ 
fective temperature, fcsT = 4.0 x 10“ 21 kg • m 2 • s~ 2 ; and 
tension of the molecular motors acting on the crossing fil¬ 
aments, T m = 2.4 x 10 -12 kg ■ m • s -2 . These parameter 
values are in a physiologically plausible range and are 


III. COLLECTIVE MOTIONS OF AFCS 

Several types of collective behaviors were observed 
depending on the shape angle 9 a and density p = 
N/L 2 p,m~ 2 . In Fig. [2j examples of observed dynam¬ 
ics are shown in snapshots at the indicated parameters 
values. In the figures, each AFC is colored according 
to its direction ipi (i = 1,2, ,N) defined as the bi¬ 

secting direction of the two filaments [Fig.[5][a) and (b)]. 
As explained below, we measured several quantities that 
characterize the AFCs’ collective behaviors and classi¬ 
fied the dynamics into four phases: globally ferromag¬ 
netic [Fig. 0[c)], moving smectic [Fig. [2][d)], disorder 


(a) Direction of an AFC 




(c) Globally ferromagnetic 
($ a = 0°, p = 5(pm" 2 )) 


(d) Moving smectic 
($ a = 30°, p = 5) 




(e) Disorder 

(@ a = 90> = 5) 


(f) Locally nematic 
(@ a =150> = 5) 




FIG. 2. Snapshots of simulation results, (a) Direction of an 
AFC, ipi, is defined as the direction bisecting two filaments 
for each AFC. (b) Color code for the direction of AFCs. (c)- 
(f) Snapshots of the simulation results at indicated parameter 
values. Filaments are colored according to their directions. 
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[Fig. [2je)] , and locally nematic [Fig. |2[f)] patterns. 


A. Characterization of AFC dynamics 


The most significant quantity for characterization of 
AFC dynamics is the orientation order parameters. We 
measured ferromagnetic (polar) and nematic (apolar) pa¬ 
rameters, R\ and f? 2 , defined by R\ = -k 


EZi ex P (iV’i) 


_ 1 


EZi ex P (2iV>i) 


, respectively. 


and f ?2 = jj- 

In some parameter regions, the system does not exhibit 
global order but rather shows strong directional align¬ 
ment among neighboring AFCs [Fig. |5[f)]. To detect 
such local structures, we divided the system into regular 
10 x 10 lattices with filament length size (t x £ [/jm 2 ]), 
measured the ferromagnetic and nematic order param¬ 
eters in each region, and then averaged their values. 
These local order parameters R l ° c and R l 2 c are given 

as R[ oc = iSoE}™ EtiexP^i) and R l 2 oc = 

,N- 


TM Ej=i Wj Ei=i ex P ( 2i ^) > respectively. Here, Nj in¬ 
dicates the number of AFCs whose anchor is in the j-th 
lattice (j = 1,2, ••• ,100). These quantities are shown 
on a color scale in Fig. [3[a)-(d). 

As we will see below, the system shows another spatial 
regularity in the moving smectic phase where the density 
bands of the aligned AFCs are observed [Fig. |2][c)]. To 



FIG. 3. (a)-(d) Color maps of R 1} R 2 , R[ oc , and R l 2 oc on 0* 
and N. (e) Definition of ips. (f) Color map of S = (sin 2 ips). 
The parameter values were calculated by the time average 
between 900 s < t < 1000 s. 


characterize the density bands in the AFCs’ positions, 
an additional order parameter, S , was introduced as the 
average of sin 2 ips, where ips is the angle defined for two 
contacting AFCs, as indicated in Fig. 0(e). If two AFCs 
are aligned in parallel and are positioned side by side, 
i.e., the line connecting their anchors is perpendicular to 
their directions (ips ~ ±90°), S becomes high. The color 
map of S is shown in Fig. [3jf ). 

To confirm that our analyses are independent of the 
system size, we conducted numerical simulations with 
different system sizes. The results are summarized in 
Appendix [A] 


B. Dynamics of AFCs 

(a) Globally ferromagnetic order. For a small angle 
between the two filaments of an AFC (0° < 0 a < 80°), R\ 
becomes nonzero as the AFC’s density increases, which 
indicates a globally ferromagnetic order [Fig. [5[c) and 
Fig. |3[a)]. R l f c also becomes nonzero at the same den¬ 
sity; hence, once the order emerges, it grows globally, 
indicating a phase transition [Fig. [3fc)] . Note that at 
6 a = 0°, an AFC is equivalent to a single polar filament, 
for which ferromagnetic transition was reported in earlier 
studies [H, HH- 

When the concentration of AFCs is high, directional 
order develops up to the system size, and defects such as 
vortices and asters are not observed. It is possible that 
these defects would be observed in the lower density re¬ 
gion near the transition point; however, it is practically 
difficult to detect these structures at the resolution re¬ 
quired for this type of simulation. Therefore, we did not 
consider vortices and asters further in the present study. 

(b) Moving smectic structure. Inside the region of 
the ferromagnetic order phase, lamellar structures that 
consist of AFCs facing the same direction are observed 
(Fig. (2[d)). The lamellae move through the system in 
one direction [Fig. 0]. We call this dynamic pattern a 
‘moving smectic’ by analogy with the physics of liquid 
crystals. Although moving band structures of density 
waves were reported in some actively propelling systems 
lii, the smectic pattern found here is distinct, 
because the dominant propelling direction of each parti¬ 
cle is nearly perpendicular to the band (See Sec. HVI) . As 
shown in Fig. [3[f), the order parameter S distinguishes 
this pattern from the globally ferromagnetic order where 
the positions of the AFCs exhibit no regularity. 

(c) Disorder. When the density of the AFCs is small, 
diffusion is dominant since the AFCs are sparse and thus 
interactions are rare. There is neither orientational nor 
spatial order. Even when the density of the AFCs is high, 
in the region 80° < 9 a < 120°, no sign of order in orien¬ 
tation, space, and time is detected [Fig. [5](e), Fig. Elf)]. 

(d) Locally nematic structure. For larger 9 a up to 180°, 
R l 2 c becomes nonzero, while R 2 remains zero regardless 
of the AFCs’ density [Fig. [3[b) and (d)]. The snapshot 
in Fig. [2[f) clearly indicates that neighboring AFCs are 
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aligned and bundled, whereas orientational order is not 
maintained at a scale larger than the length of a few 
filaments. 


C. Phase diagram 

A phase diagram of the observed dynamics is depicted 
in Fig. [5j For the high-density region (p > 3pm~ 2 ), 
disorder/ferromagnetic transition is observed with small 
6a, 0° < 6 a < 80°. A moving smectic structure ap¬ 
pears in the higher density region of 30° < 6 a < 80°. 
As 8 a increases, these two patterns vanish and the mo¬ 
tions become disorder. For large obtuse angles of 0 a , 
locally nematic orders arise in 120° < d a < 180°. The 
detailed characteristics of each pattern are summarized 
in Table [T] In the remainder of this paper, the density of 
AFCs is fixed at p = 5 /zm , where all the dynamics 
are observed by changing 9 a . The dynamics of globally 
ferromagnetic, moving smectic, disorder, and locally ne¬ 
matic patterns are analyzed by setting 9 a = 0°,30°,90°, 
and 150°, respectively (represented as stars in Fig.0. 


IV. INDIVIDUAL MOTION OF AFCS 

To further understand the collective dynamics of AFCs 
described in the preceding section, the individual mo¬ 
tions of the AFCs were studied by tracking their trajec¬ 
tories. Employing the color code indicated in Fig. [6](a), 
the trajectories of 10 AFCs are shown in each dynamic 
[Fig. (6](c)-(f), left column]. Using the trajectories, we 
measured the propelling direction of each AFC, ip v , de- 


t= 997.50 (s) t= 997.54 



t = 997.58 



t= 997.62 



FIG. 4. Time evolution of the moving smectic structure ( 8 a = 
30° and p = 5 pm~ 2 ). Two lamellae are indicated by green 
and cyan arrowheads, respectively. One lamella is traced with 
a pink line. 


termined by the positional shift of the anchor between 
two time points t and t + At. The value At is set to be 
0.01 s, during which the density pattern seldom moves 
and the AFCs more or less maintain the same contacting 
pairs [Fig. [4], Then, the relative propelling direction of 
an AFC is defined by f/v = i’v — il>i [Fig. Ob)]. The fre¬ 
quency and mean velocity of an AFC in a given direction 
if r are summarized in Fig. (5](c)-(f) (right column). 

(a) Globally ferromagnetic order. In the parameter 
region of the globally ferromagnetic order, the motion 
of an AFC is translational without significant rotational 
motion [Fig. (Bfc)]. Statistical analysis showed that they 
barely moved in the same direction as that of the AFCs. 
Most frequently, the AFCs move almost transversely to 
their direction (ip r ~ ±105°) where the mean velocity is 
fastest. The less frequently observed AFCs reversed their 
direction of movement (ip r ~ 180°) at slower velocity. 

(b) Moving smectic structure. The direction and mag¬ 
nitude of the translational motion of AFCs are similar 
to those observed in the globally ferromagnetic phase 
[Fig. Hd)]. In the moving smectic phase, where a large 
proportion of the AFCs are aligned in parallel in the lay¬ 
ered bands of the lamellar structure, the statistics indi¬ 
cated that the AFCs are moving along the density band, 
and the density band itself moves slowly in the direction 
opposite to that of the AFCs [Fig. 0]. 

These behaviors are apparently contradictory to the 
momentum conservation that must be satisfied in the 
present model, except for the noise effect. For a more 
detailed analysis, the joint distribution and mean veloc¬ 
ity on ipi and if> v were measured at a given time point and 
were compared. The joint distribution shown in Fig.[7][a) 
revealed two peaks at (if>i, ip v ) ~ (60°, —45°), (40°, 145°), 
consistent with the high frequency of i /y = — ifi 

around ±105° [Fig. (51(d)]. On the other hand, the 
magnitude of the mean velocity showed two peaks at 



(a) Globally 
ferromagnetic 

(b) Moving 
smectic 

(c) Disorder 

(d) Locally 
nematic 


o° 


30 ° 60 ° 


90 ° 

9 a 


120 ° 150 ° 180 ° 


FIG. 5. The phase diagram of the collective patterns of AFCs. 
(a) Globally ferromagnetic, (b) moving smectic, (c) disorder, 
(d) locally nematic. The precise criterion for assignment of 
phases in the numerical simulation is based on time-averaged 
order parameters as follows: globally ferromagnetic: Ri > 0.5 
and S < 0.55; moving smectic: Ri > 0.5 and S > 0.55; lo¬ 
cally nematic: Ri < 0.5 and R l f c > 0.5; disorder: otherwise. 
Stars in the figure indicate the parameters used for analyzing 
dynamics in the respective phases. The square indicates the 
parameters used in Fig. 1101 
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(V’i.V’t/) - (110°, 0°), (-40°, 70°) [Fig. mb)]. The dis¬ 
crepancy of the peak positions between frequency and 
velocity support the following scenario: although a large 
portion of the AFCs makes up the density bands that 
determine their collective behaviors, a smaller portion of 
the AFCs exhibit motion that is far from their collective 
average, which is much faster and in different directions, 
thereby compensating for the conservation law. Indeed, 
AFCs that move fast and escape from a density band 
were observed in the simulation [Fig. [2J[d)]. Such ‘divi¬ 
sion of labor’ in AFC assembly is likely responsible for 
the appearance of a moving density band in the system 
with momentum conservation. 

(c) Disorder. In the disorder phase, the motion is 
still translational, although the predominant propelling 
direction now coincides with the direction opposite to 
that of the AFCs [Fig. [6][c)] . 

(d) Locally nematic structure. The motion in the lo¬ 
cally nematic region is rotational rather than transla¬ 
tional [Fig. Elf)]. Compared to the other phases, the 
velocity of translational motion is low. Without signifi¬ 
cant translation of AFCs, the alignment process of AFCs 
occurs locally, resulting in the bundled structure of AFCs 
[Fig. Gif)]- 


V. TWO-BODY INTERACTION 

Investigating the interaction between AFCs is crucial 
for understanding the emergence of various dynamics at 
different parameter regions. In this section, we assess 
whether the completely overlapped state of two AFCs is 
sustained during the time evolution of the model. Linear 
analysis is not useful in this case due to the condition that 
interaction works only when filaments are crossing; thus, 
the state is singular. Instead, we numerically checked 
the stability of the overlapped state by measuring the fi¬ 
nal distance between the two initially overlapped AFCs. 
The distance of the two AFCs is defined by the sum of 
distances between the centers of L-L and R-R filaments 
[d = dL + dji pm, see Fig. [5](a)], by which d = 0 pm 
corresponds to the overlapped state. We measured the 
eventual distance of d pm after simulating the two-AFCs 
system with a simulation time up to t = 1000 s. By iter¬ 
ating 10,000 independent runs of the simulations starting 
from an almost overlapped but slightly perturbed state 
(d < 0.01), the distribution of the distance was evaluated 
for a given value of 9 a . The results were robust to the 
choice of the simulation time as long as t > 100 s. 

By calculating the distribution for each value of 9 ai we 
obtained the contour map of this distribution, as shown 
in Fig. [5] Eventual states are classified into six domains, 
(i)-(vi), corresponding to the configurations of the two 
AFCs depicted in Fig. (EJc). Configurations (i) and (ii) 
represent the almost overlapped and well-aligned states, 
where the distance d is smaller than 1 pm, which is equal 
to half of the filament length £/2. 

In 9 a < 10°, two types of configurations, (i) and (iii), 


are possible; however, the probability of type (iii) be¬ 
comes higher as 9 a increases. For 30° < 9 a < 80°, there 
are no domains with small d , indicating that the over¬ 
lapped state is not maintained. In particular, configu¬ 
ration (iv) dominates the probability up to 9 a < 120°. 


(a) 


(b) 


time course 


f = 997.5 (s) 997.6 



(d) Moving smectic (8 a = 30°) 




0 ° 



(e) Disorder (t9 a = 90°) 



(f) Locally nematic (0 a = 150°) 



0 ° 



0 ° 



FIG. 6. Trajectories of individual AFCs in the collective mo¬ 
tion. (a) Color code of trajectories. AFCs are colored along 
the time course from yellow to blue, and finally red. (b) The 
definition of ipi, i />„, and ip r . (c)-(f) (left) Trajectories of in¬ 
dividual AFCs in the indicated phase. Ten AFCs are colored 
according to (a); the others are colored gray, (right) The fre¬ 
quency of t p r (blue) and the average magnitude of velocity 
with respect to i]> r (orange). 



















For 8 a > 80°, an almost overlapped configuration (ii) ap¬ 
pears, and the probability becomes higher as 8 a increases, 
whereas the probability in non-aligned configurations be¬ 
comes smaller. Within this region (105° < 9 a < 150°), 
configuration (vi), in which two AFCs are loosely aligned, 
is also possible. 

Figurellfd) (red lines) shows the ratio so that the even¬ 
tual distance d is smaller than 1 pm [i.e., domains (i) and 
(ii)]. As 8 a increases from 8 a = 0°, configurations (iii) 
and (iv) become more dominant and the ratio decreases 
to almost 0 pm ( 8 a > 10°). The overlapped state is not 
stable where two AFCs move in opposite directions and 
are repulsive, which is responsible for the translational 
motions seen in the globally ferromagnetic, moving smec¬ 
tic, and disordered states. The ratio increases with larger 
values of 8 a ■ With the disappearance of configuration 
(iv) where 9 a > 120°, the overlapped state tends to be 
sustained as configuration (ii). This observation of inter¬ 
action between two AFCs is critical for the local nematic 
pattern, in which translational motion is suppressed and 
the AFCs are well aligned. 

VI. CROSSING OF AFCS AND FILAMENTS 


arise from the variety of crossing combinations for each 
filament. Measurements of the number of crossing points 
between two contacting filaments revealed that in most 
cases, except when considering globally ferromagnetic 
dynamics, the filaments dominantly crossed at a single 
point [Fig- fTOT al]. The abundance of four-point crossings 
in the globally ferromagnetic phase is easily explained by 
the small angle of 6 a . Therefore, we only consider one- 
point crossings and classify them into four types, namely 
RR, LR, RL, and LL. The frequencies of these crossing 
types are shown in Fig. HOl' bf. 

(a) Single active filament system. First, we investi¬ 
gate crossing in a single active filament system ( 8 a =0°) 
where <p s and <p c are identical. Distributions of the an¬ 
gles are dependent on the density p [Fig. |5[b)]. Despite 
the global order at p = 5 pm -2 , the angular distribution 
P{4>s) is not highest at cj) s = 0°, indicating that contact¬ 
ing filaments are not aligned in parallel but are instead 
tilted toward each other ( (f> c ^0°); the mean and stan¬ 
dard deviation of the angular difference are calculated 
from the distribution of |^ s | to be p = 66° and a = 39°, 
respectively. 

(b) Global ferromagnetic order. In a small angle of 


Having elucidated the interaction between two fila¬ 
ments, in this section we investigate how AFCs in a pop¬ 
ulation cross each other to achieve characteristic patterns 
depending on the parameters. In particular, because only 
the angle between crossing filaments determines the de¬ 
velopment of an AFC’s position in Eq. © , our main 
interest is the relative angle distributions among inter¬ 
acting AFCs. Two kinds of angles are considered below. 
The first is the angular difference between crossing fil¬ 
aments, </> s , as shown in Fig. [9|a). The second is the 
angular difference between AFCs, <f c , measured between 
the directions of contacting AFCs [Fig. ©c)]. 

Because an AFC has two filaments, complications can 


(a) freauencv (b) mean velocity turns 1 ') 



0.2 

0.15 

0.1 

0.05 

0 


FIG. 7. Statistical analyses in the moving smectic phase, (a) 
Joint distribution of ipi and if v . The peak positions are at 

— (60°,— 45°), (40°, 145°). (b) Magnitude of the 

AFCs’ velocity on ifi and ip v . The peak positions are at 

— (110°, 0°), (—40°, 70°). The majority of the AFCs 
in the population move toward ~ —45° or 145°, which 
is insufficient for momentum to be conserved. A few AFCs 
move faster toward ip v ~ 0° or 70°, which compensates for 
the momentum conservation. 



FIG. 8. Analyses for a two-body interaction, (a) The distance 
of the two AFCs is defined by d = dL + dn, where di(dij) 
is the distance between the centers of L-(R-)filaments. (b) 
Distribution of d after the simulation for each 6 a . (i)-(vi) 

denote domains with characteristic configurations, as depicted 
in (c). (d) The ratio that the final distance of two AFCs is 
smaller than 1 pm. 
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9 a corresponding to the globally ferromagnetic phase, 
the distribution of (j> c becomes different from that of <fi s 
[Fig. EJd)]. The angular distribution P(4> c ) now exhibits 
a peak at 9 C = 0°, which means that parallel alignment 
of two AFCs are enhanced more than for a single active 
filament owing to the change in molecular shape. This 
enhancement in AFC alignment explains the increase in 
Ri against 9 a in Fig. HUA) (Appendix 0 and the lower 
critical density for larger 9 a (0° < 9 a < 80°) in Fig. [5] 

(c) Moving smectic structure. The distribution of the 
angular difference between contacting filaments ((f> s ) is 
small (n = 50°) and is sharply distributed (cr = 28°). 
With higher frequency of the LR and RL types of cross¬ 
ings compared to that of the other phases [Fig. [TOT bl]. 
the distribution of (f> c shows a significantly high peak at 
4> c = 0°, indicating a strong tendency of AFCs to be 
directed in parallel, which underlies the formation and 
maintenance of moving density bands. Note that such 
parallel alignment of AFCs is not predicted from the two- 
body analysis (See Fig. [5] in Sec. 0, and it results from 
many body interactions. 

(d) Disorder. The angle between crossing filaments, 
cj) s , is most frequently found in the perpendicular direc¬ 


tion ((f> s ~ 90°). This distribution is almost identical to 
those observed for collections of single filaments with ran¬ 
dom orientation [Fig. [9jf )]• The frequency of more-than- 
two-points crossings is small [Fig. llOf aV. Consistent with 
a random orientation, the distributions of <p c and cross¬ 
ing types [Fig. ITOl' b'l] arc almost uniform. These results 
indicate that orientational order does not exist, not even 
locally. 

(e) Locally nematic structure. In the locally nematic 
phase, local alignment processes make the AFCs overlap 
and form a bundle structure. In the bundle, the align¬ 
ment is almost perfectly overlapped, as represented by a 
sharp peak in the distribution of the angular differences 
of AFCs at </> c = 0° [Fig. EH[h)]. The frequency of more- 
than-two-points crossings is negligibly small [Fig. llOf c)]. 
In the distribution of angular differences in crossing fil¬ 
aments, several peaks appear [Fig. ESh)]. These peaks 
can be explained by interactions of two AFCs, as shown 
in Fig. (Kb) and (c). The peaks at <p s = 0° are due 
to LL/RR crossings corresponding to configuration (ii), 
whereas the peaks at </> s = ±9 a are due to LR/RL cross¬ 
ings corresponding to configuration (iii). The other two 
peaks in the neighbors of (f> s = 0° are produced by con¬ 
figuration (vi). 


(a) (b) Single active filament (8a = 0°) 




(d) Globally ferromagnetic (8a = 5°) (e) Moving smectic (8a = 30°) 



(f) Disorder (8a = 90°) (g) Locally nematic (da = 150°) 



FIG. 9. Statistics describing how the AFCs and filaments 
cross each other, (a) Definition of <j> s , the angular difference 
between crossing filaments, (b) Distribution of (f> s for single 
filaments on p = 1, 2.5, 5 /im -2 . (c) Definition of <j) c , the an¬ 
gular difference between crossing AFCs. (d)-(g) Distributions 
of 4> s and </> c are denoted by green and red lines, respectively. 


VII. DISCUSSION 


In the majority of previous studies, objects with 
simpler shapes were assessed, such as material points, 
spheres, rods, or ellipsoids, whether stiff or flexible 
[2ll . [23j . Shapes that are more complex, however, can 
play important roles in biological systems. For example, 
the self-organization of microtubules ~ 10 nm beneath 
the apical membrane of multi-ciliated cells is associated 
with a regular array of cilia [Hi, [39j. In the present study, 
with the aim of understanding the dynamics of the molec¬ 
ular complexes of active filaments, we performed two- 
dimensional simulations of interacting AFCs defined as 
connected pairs of biofilaments, whose shapes were char¬ 
acterized by the angle between the two filaments. 

The AFCs exhibited several types of collective patterns 
with global and local orders. These patterns can be un- 




(»a = 5°) (fl a = 30°) (» a = 90°) (6> a = 150“1 


Globally Moving Disorder Locally 
Ferromagnetic Smectic nematic 

(0a = 5°) (0 a = 30°) (0 a = 90°) (0 a = 150°) 


FIG. 10. Frequencies of the crossing types of AFCs. (a) Fre¬ 
quencies of the number of crossing points between two con¬ 
tacting AFCs. (b) Frequencies of crossing types. RL-crossing 
AFCs are shown in the figure. 
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TABLE I. Summary of emergent patterns (the light-anchor model). 


Patterns 

Globally 

ferromagnetic 

Moving 

smectic 

Disorder 

Locally 

nematic 

9a 

0° - 80° 

30° - 80° 

80° - 120° 

120° - 180° 

Density transition 

O 

O 

X 

X 

Characteristic 

movement 

Translational 
perpendicularly 
to AFC 

Translational 
perpendicularly 
to AFC 

Translational 
parallelly 
to AFC 

Rotational 

Alignment 
b/w crossing AFC’s 

tilted 

parallel 

None 

Strong 

2-times symmetry 

Alignment 

propagation 

Isotropic 

Perpendicular 
to AFC 

None 

None 


derstood by tracing the trajectories of AFCs and inves¬ 
tigating the alignment between crossing pairs of AFCs, 
as described in Sec. |TV] [VI] Analysis of the two-body 
interaction in Sec. revealed that the overlapped state 
of AFCs was maintained where 9 a was large, which was 
responsible for the appearance of the local nematic or¬ 
der. In contrast, instability of the overlapped state in¬ 
dicates repulsion between the two AFCs, which is re¬ 
quired for translational motion. Translational motion is 
necessary for establishing the global order, because the 
phase transition of continuous order parameters in a two- 
dimensional system is possible only when a long-range 
interaction is present j4g|. As discussed in Sec. IIVI and 
Sec. IVI1 ferromagnetic orders can actually grow over a 
range that is much longer than the filament length if the 
AFCs move translationally and the alignment process can 
propagate. 

Among the observed patterns, the moving smectic pat¬ 
tern is of particular interest. In contrast to earlier stud¬ 
ies that reported a moving density band [H, Ho|, H2 ; ' n 
our system, the total momentum is conserved when the 
noise is eliminated. The moving density band emerges 
even in the no-noise limit condition, which apparently 
contradicts with momentum conservation. As described 
in Sec. IIV1 a small population of AFCs that move fast in 
an atypical direction can compensate for the momentum, 
which leads to the appearance of a moving smectic. In 
the present system, a complex molecular shape is respon¬ 
sible for this mechanism to arise. 

The present study ignored some potentially important 
interactions for modeling experimental situations. For 
example, it is possible that the length of filaments and 
the shape of AFCs are not uniform and instead change 
during the time evolution. In particular, active regula¬ 
tion of filament length can be important for establishing 
ordered patterns, as exemplified by microtubule bundles 
in plant cells. Furthermore, excluding the volume ef¬ 
fect can result in alignment of the filaments and appear¬ 
ance of nematic order |2 tI . All of these processes can 


contribute to the organization of biofilament networks in 
actual systems. In further research, the present analy¬ 
ses for AFCs could be easily extended, and this would 
likely provide useful insights for understanding collective 
patterns in more complex scenarios, including the afore¬ 
mentioned processes. Employing continuum description 
[4ll I47l - l49| will clarify the mechanism by which these 
patterns emerge as well as their stabilities. The rapid 
developments in bio-nano engineering will also enable a 
comparison between experimental and simulated systems 

M. 
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Appendix A: System size dependence 

To check the system size dependence of the dynamics 
observed in the main text, we conducted numerical sim¬ 
ulations with different system sizes. We fixed the density 
as p = 5 and changed the system size as L = 10,20, 

and 40 /im. The results are shown in Fig. [TT] (a)-(c). R\ 
and I ?2 on acute 0 a are independent of the system size, 
and finite size scaling analysis supports that the emer¬ 
gence of ferromagnetic order is a genuine phase transi¬ 
tion, consistent with observation in rod-like filament dy¬ 
namics [30]. 

On the other hand, R\ and R 2 on obtuse 9 a decrease 
as the system size increases, suggesting that the globally 
nematic pattern is at most a quasi-long range order akin 
to an XY-model [2(| |22|. S is independent of the system 
size. 
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